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Abstract 


We  demonstrate  a  new  generalized  least-squares  fitting  method  which  can  be  used  to  estimate 
the  slope  of  the  best-fitting  straight  line  that  results  when  two  separate  data  sets  which  are 
expected  to  be  linearly  correlated  are  subject  to  different  uncertainties  in  their  measurements. 
The  algorithm  determines  not  only  the  optimum  slope,  but  also  produces  estimates  of  the 
intrinsic  errors  associated  with  each  data  set.  It  requires  almost  no  initial  information  about 
the  errors  in  each  data  set. 


Executive  Summary 


One  of  the  problems  that  frequently  confronts  an  experimentalist  is  that  of  fitting  a  straight 
line  to  data,  that  is,  the  problem  of  determining  the  functional  relationship  between  two 
variables.  A  least-squares  calculation  is  the  most  common  approach  to  this  problem  and,  in 
its  simplest  form,  one  variable  is  subject  to  error  while  the  other  is  error- free.  This  procedure 
is  described  in  almost  all  elementary  books  on  statistics. 

A  slightly  more  complex  case  is  that  both  variables  have  errors,  which  are  well  known. 
Considerable  work  has  been  done  on  fitting  procedures  which  take  into  account  errors  in 
both  variables.  These  fitting  procedures  generally  work  very  well  for  the  tasks  for  which 
they  are  assigned.  However,  there  are  times  when  two  data  sets  are  compared  in  which  the 
intrinsic  measurement  errors  of  one  or  both  variables  are  unknown.  In  this  case,  there  is 
no  algorithm  available  to  determine  the  best-fit  line.  The  main  objective  of  this  paper  is  to 
demonstrate  a  new  procedure  for  this  objective. 

In  this  paper,  we  demonstrate  a  new  generalized  least-squares  fitting  method  which  can  be 
used  to  estimate  the  slope  of  the  best-fitting  straight  line  that  results  when  two  separate 
data  sets  which  are  expected  to  be  linearly  correlated  are  subject  to  different  uncertainties 
in  their  measurements.  The  algorithm  determines  not  only  the  optimum  slope,  but  also 
produces  estimates  of  the  intrinsic  errors  associated  with  each  data  set.  It  requires  almost 
no  initial  information  about  the  errors  in  each  data  set.  The  algorithm  works  best  if  the 
user  knows  which  of  the  two  variables  has  the  smallest  intrinsic  error,  but  is  not  limited  to 
this  condition.  We  hope  that  this  process  will  find  application  in  a  wide  range  of  situations, 
e.g.,  ranging  from  radar  and  satellite  measurements  to  any  physical  measurements. 
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1  Introduction 


Least-squares  techniques  for  examining  best-fit  lines  to  correlated  data  sets  are  well  estab¬ 
lished,  but  in  general  are  restricted  to  cases  in  which  the  errors  in  the  two  data  sets  being 
compared  are  known.  The  simplest  case  is  that  in  which  there  is  no  error  in  one  variable 
(usually  plotted  as  the  abscissa),  and  the  error  in  the  other  variable  is  known.  This  proce¬ 
dure  is  described  in  almost  all  elementary  books  on  statistics  [e.g.,  Young,  1962;  Bevington, 
1969;  Daniel  and  Wood,  1971;  Brandt,  1976;  Taylor,  1982]. 

A  slightly  more  complex  case  is  that  in  which  both  variables  have  errors,  but  both  errors  are 
well  known.  Examples  of  such  process  have  been  shown  by  (amongst  others)  York  [1966], 
Barker  [1974],  Orear  [1982],  Lybanon  [1984],  Miller  and  Dunn  [1988],  Reed  [1989,  1992],  and 
J olivette  [1993].  These  fitting  procedures  generally  work  very  well  for  the  tasks  for  which 
they  are  assigned. 

However,  there  are  times  when  two  data  sets  are  compared  in  which  the  intrinsic  measure¬ 
ment  errors  of  one  or  both  variables  are  unknown.  In  this  case,  there  is  no  algorithm  available 
to  determine  the  best-fit  line.  The  purpose  of  this  paper  is  to  demonstrate  a  new  procedure 
for  this  objective. 

This  paper  is  broken  up  as  follows.  We  begin  by  simulating  some  data  using  computer 
techniques,  in  which  the  intrinsic  errors  are  pre-specified.  We  then  demonstrate  the  effects  of 
employing  existing  methods  of  data  analysis  to  these  data  sets,  and  point  out  the  limitations 
of  these  procedures.  We  then  demonstrate  how  these  standard  processes  lead  to  our  new 
procedure,  and  develop  our  new  algorithm.  Finally,  we  use  new  artificially  generated  data 
to  demonstrate  the  application  of  our  new  method. 


2  Standard  analysis  procedures 


As  an  initial  study,  we  used  computer-generated  data  to  re-visit  standard  least-squares  fitting 
procedures.  Our  model  works  as  follows.  We  generate  two  data  sets,  one  representing  the 
abscissa  (xt),  and  the  other  the  ordinate  (yt).  These  are  generated  by  randomly  selecting 
values  Xj  from  a  Gaussian  distribution  with  a  pre-specified  standard  deviation.  In  our  first 
case,  we  will  demonstrate  the  situation  for  a  standard  deviation  of  Y,x  =  40  and  Yy  =  40. 
We  then  calculated  y,-  values  using  the  expression 
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Vi  =  9o  Xi  +  Cl 


(1) 


For  our  initial  calculations,  we  have  chosen  g0  =  1,  which  produces  a  standard  deviation 
in  the  ?/,-  values  (Sy)  equal  to  Sx.  Notice  that  these  2  values  are  not  errors,  but  represent 
the  spread  in  our  raw  data  points.  Other  options  are  considered  later.  This  process  then 
gave  a  straight  line  graph.  Our  next  step  was  to  partly  randomize  the  points.  For  each 
point,  we  added  a  random  number  to  the  x>  value,  and  a  different  random  number  to  the 
yt  value,  where  each  randomly  generated  number  was  derived  from  a  Gaussian  distribution 
with  standard  deviations  ax  and  ay  respectively.  We  generated  about  300  points  per  such 
realization  in  the  first  instance. 

Our  next  step  was  to  perform  various  types  of  least-squares  fitting  procedures.  We  removed 
the  mean  x  value  from  the  x{  values,  and  the  mean  y  value  from  the  t /,•  values,  and  then 
applied  standard  fitting  procedures  in  two  different  ways.  Firstly,  we  fitted  the  data  assuming 
that  all  the  X{  values  have  no  error.  This  gave  an  equation  which  we  denote  as 


Vi  =  9x  Xi  +  c2  (2) 

Because  the  data  were  normally  distributed,  as  were  the  errors,  we  found  that  c2  was  generally 
close  to  zero.  We  will  therefore  not  further  investigate  this  value,  and  will  concentrate  on  the 
parameter  gx.  For  real  data,  we  assume  that  the  user  will  have  already  removed  the  means 
before  applying  our  method. 


Next  we  fitted  the  same  data  assuming  that  the  yi  values  have  no  error,  and  obtained  lines 
of  the  type 


Xi  =  (1  /gy)yi  +  c3 


(3) 


or 


Vi  -  gy  Xi  +  c4 


(4) 


As  a  final  check,  we  also  fitted  our  data  using  the  algorithm  described  by  York  [1966],  in 
which  we  used  our  known  values  for  ax  and  oy  to  estimate  the  mean  slope.  We  will  denote  this 
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Figure  1:  Examples  of  the  slopes  obtained  from  three  different  methods  applied  to  some 
sample  data.  The  symbols  gx  and  gy  refer  to  the  slopes  obtained  from  the  regressions  y  on  x 
and  x  on  y  respectively,  while  gxy  refers  to  the  slope  obtained  from  the  algorithm  described 
by  York  [1966].  crx  and  <ry  are  measurement  errors  in  x  and  y  respectively. 

estimate  as  gxy.  This  latter  step  is  not  required,  but  was  simply  done  as  a  consistency  check. 
We  then  repeated  this  process  for  each  ( <jx ,  ay)  pair  for  many  different  realizations  (typically 
300).  Figure  1  shows  examples  of  our  fitting  attempts  for  cases  with  (crx  =  5,  cry  =  10)  and 
(<rx  =  5,  <7y  =  20). 

As  noted,  we  have  repeated  such  fits  for  many  hundreds  of  different  realizations,  and  for 
different  combinations  of  ax  and  ay,  and  different  numbers  of  points  per  sample.  The  results 
of  the  mean  slopes  gx,  gy,  and  gxy,  and  associated  standard  deviations  for  the  mean  (the 
latter  being  denoted  as  A gx,  Agy,  and  A gxy)  are  shown  as  contour  plots  in  Figure  2.  It  is 
clear  that  the  York  [1966]  algorithm  has  worked  well  (bottom  panel),  with  all  slope  estimates 
being  close  to  unity  (the  original  slope).  We  will  therefore  not  further  discuss  the  results  of 
the  York  algorithm. 

However,  the  distributions  of  gx  and  gy  are  striking  in  the  way  that  the  contours  are  aligned 
almost  parallel  to  the  axes  (see  Figure  2).  It  is  this  rather  special  alignment  which  permits 
us  to  develop  the  algorithm  described  in  this  paper.  The  arrangement  of  these  lines  shows 
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Figure  2:  Slopes  obtained  from  three  different  methods  as  a  function  ax  and  cry.  gx  and  gy 
refer  to  the  slopes  obtained  from  the  regressions  y  on  x  and  x  on  y  respectively.  gxy  refers 
to  the  slope  obtained  from  the  algorithm  described  by  York  [1966].  <jx  and  ay  are  intrinsic 
noise  and  measurement  errors  in  x  and  y  respectively.  A gx,  Agy,  and  Agxy  are  the  standard 
deviations  for  the  means  (standard  error)  of  these  various  slopes. 


that  the  estimate  of  gx,  although  erroneous,  is  pretty  much  independent  of  ay.  Likewise,  gy 
depends  very  weakly  on  ax.  Thus  we  may  write 


9x  ~  9x  M  ;  9y~gy  ( cry )  (5) 

There  are  slight  deviations  from  this  law  for  small  ax  and  large  ay  in  the  top  graph,  but 
otherwise  these  dependencies  are  fairly  true.  Because  of  this  striking  character,  we  then 
decided  to  develop  an  analytical  expression  for  gx  as  a  function  of  <jx.  Our  results  are  shown 
in  Figure  3,  where  we  plot  gx  as  a  function  of  ax  for  selected  ax.  We  have  actually  taken  the 
case  ay  =  0,  but  as  noted  gx  depends  only  very  weakly  on  <jy.  We  have  also  re-plotted  our 
ordinate  and  abscissa  as  ax/T,x  and  gx/g0.  We  have  done  this  because  we  expect  that  these 
scalings  should  apply,  and  we  have  confirmed  this  by  repeating  our  analyses  using  different 
values  of  gQ  and  £*.  Henceforth  we  will  generally  work  with  such  normalized  values. 

At  this  point,  we  need  to  indicate  an  important  assumption  concerning  the  sign  of  the  slope. 
For  our  preceding  examples,  we  used  cases  for  which  the  best  fit  lines  between  data  sets  had 
a  positive  slope.  This  may  seem  restrictive,  but  in  fact  it  is  not,  as  will  now  be  seen.  In 
any  realistic  situation  in  which  the  data  show  reasonable  correlation,  gx  and  gy  will  have  the 
same  sign,  so  it  is  easily  determined  what  the  sign  of  the  true  best  fit-slope,  gQ ,  should  be.  If 
it  is  negative,  we  can  simply  reverse  the  signs  of  all  of  the  ?/t  values,  and  the  best  fit  line  will 
then  have  positive  slope.  We  may  then  apply  the  algorithm,  and  upon  its  completion  (and 
determination  of  the  best  fit  slope  g0 ),  we  may  then  get  our  true  best-fit  slope  by  reversing 
the  sign  of  g0.  Thus  our  restriction  to  positively  correlated  data  is  not  a  limitation  of  the 
technique,  and  the  rest  of  the  discussion  will  proceed  assuming  such  positively  correlated 
data. 
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Figure  3:  Variation  of  crx/T,x  as  a  function  of  gx/go ■  The  solid  line  with  asterisk  symbols 
is  the  empirical  formula  given  by  equation  (6),  and  the  diamond  symbols  are  the  values 
calculated  from  direct  application  of  the  numerical  model  described  in  Section  2. 


We  may  now  return  to  our  developmental  discussion.  Following  generation  of  Figure  3,  we 
determined  an  analytical  function  which  fitted  these  points.  The  curve  is  continuous  and 
continuously  differentiable,  and  is  plotted  as  the  solid  line  in  the  figure.  The  expression  we 
have  used  is  complicated,  but  quite  accurate,  and  is  given  by 
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(6) 


A  more  succinct  formula  may  be  found  in  due  course,  but  for  now  our  purpose  is  just  to 
demonstrate  our  method.  Furthermore,  we  also  recognize  that  the  relation  between  g0/gy 
and  ay  obeys  the  same  expression  (although  note  that  we  use  g0lgy  here  compared  to  gx/g0 
above). 


2.15  Gy 
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i.e.,  as  above  with  gxjg0  replaced  by  ( gy/g0 )  1- 

This  now  means  that  if  we  have  a  set  of  points  (xj,z/,),  and  we  know  the  intrinsic  error  in 
one  of  these  variables,  we  can  determine  uniquely  the  “true”  best  fit  line  gQ  and  also  the 
intrinsic  error  in  the  other  variable.  For  example,  suppose  the  user  knows  gx.  Then,  one 
can  calculate  the  experimental  standard  deviations  gx  and  qy  for  the  two  data  sets.  These 
standard  deviations  include  both  the  natural  spread  of  the  data  and  the  intrinsic  errors.  Then 
we  may  determine  an  estimator  for  £x,  which  we  will  denote  as  £*,  through  the  relation 
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(8) 


=  \Ai?  -  ^ 

Next,  the  user  applies  standard  fitting  procedures,  first  assuming  crx  —  0  to  get  gx ,  and  then 
o-j,  =  0  to  get  gy.  One  may  use  Figure  3  or  equation  (6)  to  determine  an  estimator  for  gx/g0 
(which  we  will  denote  as  g~xo)  from  the  knowledge  of  ax/£x.  Thus,  knowing  gx,  and  gx/g0, 
one  can  readily  determine  g0,  our  estimate  of  g0.  Furthermore,  now  that  g0  is  known,  the 
user  can  determine  an  estimator  for  gylg0  (which  we  will  denote  as  g~yo).  Further  application 
of  equation  (7)  then  permits  determination  of  an  estimator  for  ay/ Ey.  This  may  in  turn  be 
used  to  determine  £y  through  the  relation 
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(9) 


Since  <;y  d"  an  estimator  for  ay  can  be  uniquely  determined  as  well.  Hence  given 

ax,  the  user  can  determine  the  slope  of  the  best-fit  line,  g0,  and  also  determine  dy.  This  is 
useful,  but  it  is  not  the  full  story. 

The  above  description  assumes  that  the  user  knows  the  intrinsic  error  <jx.  But  what  if  this 
is  not  known  ?  Can  we  determine  it  ? 


Suppose  we  “guess”  ax ,  and  carry  out  the  above  process.  Then  we  can  obtain  a„,  XL  E 
and  dy  for  this  proposed  value  of  ax.  However,  we  do  not  know  if  we  have  chosen  the  correct 
value  for  ax.  This  dilemma  can  be  solved,  because  there  is  one  extra  piece  of  information 
available  which  has  not  yet  been  used.  The  means  have  been  removed,  so  that  if  there  were 
no  random  errors,  we  would  expect  Zy/Zx  should  equal  g0.  Hence  we  may  say  that,  to  quite 
high  accuracy, 


Ey/Ea,  ~  g0  (10) 

It  is  essential  that  this  condition  is  closely  satisfied  in  order  for  our  selection  of  ax  to  be 
considered  as  acceptable.  We  may  thus  use  this  as  our  final  diagnostic  test.  If  we  try 
different  values  of  crx,  and  repeat  the  above  process,  we  can  examine  the  ratio  (ty/£x)lg0, 
and  choose  the  value  of  ax  which  allows  this  parameter  to  be  closest  to  1.  Alternatively,  we 
may  find  the  value  where  the  parameter  E  =  [((E^/E x)/g0)  -  1]  is  closest  to  zero. 
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A. 

Remove  the  means  from  both  data  sets. 

B. 

Do  standard  least-squares  fit  to  obtain  gx  and  gy. 

C. 

If  gx  and  gy  are  negative,  reverse  the  signs  on  all  y,  values, 
and  the  signs  of  gx  and  gy. 

D. 

Find  experimental  standard  deviations  of  the  data  sets,  qx  and  q 
(Note  these  include  both  natural  spread  and  intrinsic  errors). 

E. 

Determine  appropriate  steps  for  ax,  (let  these  step  sizes  be  denoted  as  Acrx). 

F. 

Step  <j x  from  zero,  upward  in  steps  of  A ax. 

G. 

For  each  value  of  crX)  do  the  following. 

Gl. 

Find  E x  =  yJ$-  a2x. 

G2. 

Find  ax/Tix,  and  then  use  equation  (6)  to  determine  an  estimate 
for  gx/g0  (call  it  g~xo). 

G3. 

From  g~xo,  and  from  known  value  of  gx,  determine  g0  (an  estimate  of  g0). 

G4. 

From  g0 ,  and  the  previously  determined  value  of  gy,  determine  g~yo  =  gy/g0. 

G5. 

Determine  Ry  =  ay/  Ey  from  equation  (7). 

G6. 

Determine  Ey  from  Ey  =  <;y/\fl  +  Ry2 . 

G7. 

Determine  ay  =  J \2  —  Sy2 . 

G8. 

V  ^  ^ 

Now  calculate  E  =  ((E y/T,x)/g0)  —  1. 

G9. 

Go  to  Gl,  and  repeat  the  loop  Gl  to  G8  and  search  for  the  point 
where  E  is  closest  to  zero. 

H. 

If  a  sign  reversal  was  required  in  step  C,  reverse  the  signs  of  all 
the  \ji  values  again,  and  reverse  the  sign  of  the  resulting  slope  g0. 

I. 

Add  the  means  back  to  the  resulting  equation. 

Note:  We  recommend  choosing  the  {xi}  and  {y,}  values  so  that  crx  <  ay. 

Table  1:  Procedure  for  unique  determination  of  g0 ,  crx,  and  cry. 

Thus  repeated  application  of  successive  values  of  ax,  coupled  with  the  above  test  for  E  =  0, 
permit  us  to  determine  the  best  estimator  for  ax.  Hence  this  iterative  approach  allows  us 
to  uniquely  determine  estimates  for  g0,  crx,  and  ay.  No  prior  knowledge  about  any  of  these 
parameters  is  required,  although  we  will  show  shortly  that  from  a  pragmatic  perspective  it 
is  wisest  to  choose  the  variables  {xj}  and  {y,}  such  that  ax  <  ay.  Our  basic  procedure  is 
summarized  in  Table  1. 
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g0  Ey/Sx  2 

0.50  38.62  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.984  16  9 

1.00  38.61  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0  984  16  6 

1.50  38.59  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.984  16.2 

2.00  38.57  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.985  15  6 

2.50  38.54  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.985  14.8 

3.00  38.50  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.986  13.9 

3.50  38.46  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0  987  12  8 

4.00  38.41  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0  989  11  5 

4.50  38.36  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.990  10  0 

5.00  38.30  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0  992  8  4 

5.50  38.23  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0  994  6  6 

6.00  38.15  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0.996  4.6 

6.50  38.07  38.62  16.43  37.98  41.38  0.96  1.19  1.000  0  998  2  4 

7.00  37.98  38.62  16.39  38.00  41.38  0.96  1.19  1.001  1.000  1.4 

7.50  37.89  38.62  16.20  38.08  41.38  0.96  1.19  1.005  1.005  0  3 

8.00  37.78  38.62  15.99  38.17  41.38  0.96  1.19  1.009  1.010  0.0 

8.50  37.67  38.62  15.76  38.26  41.38  0.96  1.19  1.014  1.016  -0.8 

9.00  37.56  38.62  15.51  38.36  41.38  0.96  1.19  1.019  1  021  -1  9 

9.50  37.43  38.62  15.22  38.48  41.38  0.96  1.19  1.025  1  028  -2  8 

10.00  37.30  38.62  14.87  38.62  41.38  0.96  1.19  1.032  1.035  -3  1 

10.50  37.17  38.62  14.43  38.78  41.38  0.96  1.19  1.040  1  044  -3  0 

11.00  37.02  38.62  13.87  38.99  41.38  0.96  1.19  1.051  1.053  -2  0 

11.50  36.87  38.62  13.26  39.20  41.38  0.96  1.19  1.063  1.063  -2.6 

12.00  36.71  38.62  12.70  39.38  41.38  0.96  1.19  1.074  1.073  -2.8 

12.50  36.54  38.62  12.20  39.54  41.38  0.96  1.19  1.084  1  082  -3  1 

13.00  36.37  38.62  11.73  39.68  41.38  0.96  1.19  1.095  1.091  -3.3 

13.50  36.18  38.62  11.21  39.83  41.38  0.96  1.19  1.105  1.101  -3.6 

14.00  35.99  38.62  10.53  40.02  41.38  0.96  1.19  1.117  1.112  -4.2 

14.50  35.80  38.62  9.59  40.25  41.38  0.96  1.19  1.128  1.125  -4.3 

15.00  35.59  38.62  8.36  40.53  41.38  0.96  1.19  1.141  1.139  -4.6 


Table  2:  Sequence  of  calculations  for  the  parameter  S  when  input  values  are  T,x  =  40  ax  =  8 
Ey  =  40,  and  ay  =  16. 


3  Tests  of  our  models 


Having  established  our  numerical  iterative  model  for  g0,  ax,  and  <fy,  it  is  of  course  necessary 
to  test  its  accuracy.  For  this,  we  again  turn  to  simulation. 


We  repeated  the  process  described  in  section  2  to  generate  our  original  {.t,  }  and  {?/,}  values, 
but  this  time  for  a  wide  variety  of  values  for  g0,  crx,  and  ay.  Then,  for  each  data  set,  we 
applied  the  algorithm  described  in  Table  1.  Table  2  shows  a  sequence  of  calculations  for  ax 
in  steps  of  0.5  (un-normalized  units).  The  important  parameter  to  observe  is  E,  as  written 
in  the  last  column  of  the  table.  In  this  particular  example,  we  see  that  E  passes  through  zero 
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Ej. 

Sr 

<Ty 

Ey 

*>y 

9x 

9y 

9o 

Sy/S* 

0.50 

40.99 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.944 

28.8 

1.00 

40.98 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.945 

28.6 

1.50 

40.97 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.945 

28.2 

2.00 

40.95 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.946 

27.6 

2.50 

40.92 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.948 

27.0 

3.00 

40.88 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.949 

26.1 

3.50 

40.84 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.951 

25.1 

4.00 

40.80 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.953 

24.0 

4.50 

40.75 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.956 

22.7 

5.00 

40.69 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.959 

21.2 

5.50 

40.62 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.962 

19.6 

6.00 

40.55 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.965 

17.8 

6.50 

40.48 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.969 

15.9 

7.00 

40.39 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.973 

13.8 

7.50 

40.30 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.977 

11.5 

8.00 

40.21 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.982 

9.1 

8.50 

40.10 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.987 

6.5 

9.00 

39.99 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.993 

3.8 

9.50 

39.88 

40.99 

21.39 

79.69 

82.51 

1.89 

2.14 

2.000 

1.998 

0.9 

10.00 

39.76 

40.99 

20.72 

79.87 

82.51 

1.89 

2.14 

2.009 

2.009 

0.1 

10.50 

39.63 

40.99 

19.63 

80.14 

82.51 

1.89 

2.14 

2.021 

2.022 

-0.7 

11.00 

39.49 

40.99 

18.11 

80.50 

82.51 

1.89 

2.14 

2.036 

2.038 

-1.2 

11.50 

39.35 

40.99 

15.97 

80.95 

82.51 

1.89 

2.14 

2.054 

2.057 

-1.4 

12.00 

39.20 

40.99 

13.36 

81.42 

82.51 

1.89 

2.14 

2.076 

2.077 

-1.8 

12.50 

39.04 

40.99 

10.75 

81.81 

82.51 

1.89 

2.14 

2.097 

2.095 

-2.6 

13.00 

38.88 

40.99 

8.31 

82.09 

82.51 

1.89 

2.14 

2.117 

2.112 

-2.8 

13.50 

38.71 

40.99 

5.93 

82.30 

82.51 

1.89 

2.14 

2.136 

2.126 

-4.6 

14.00 

38.53 

40.99 

3.50 

82.44 

82.51 

1.89 

2.14 

2.155 

2.140 

-7.4 

14.50 

38.34 

40.99 

0.89 

82.51 

82.51 

1.89 

2.14 

2.175 

2.152 

-10.9 

15.00 

38.15 

40.99 

1.96 

82.49 

82.51 

1.89 

2.14 

2.196 

2.162 

-15.7 

Table  3:  Sequence  of  calculations  for  the  parameter  S  when  input  values  are  T,x  =  40, 
ax  =  10,  =  80,  and  cry  =  20. 


at  ox  —  8.0  (highlighted  in  italics).  We  see  that  our  program  estimates  cry  to  be  15.99  and 
g0  to  be  1.009.  Our  estimates  are  very  close  to  our  original  input  values  ax  —  8.0,  ay  =  16.0, 
and  gQ  =  1.0.  The  fact  that  our  best  estimate  for  g0  is  about  1%  higher  than  the  true  value 
is  a  concern,  but  not  a  major  one.  There  are  clearly  cases  which  give  a  best  fit  g0  which  is 
better  than  our  final  estimate,  especially  at  low  ax  and  high  ay,  but  these  have  erroneous 
estimates  of  ax  and  oy.  Our  final  estimate  gives  the  best  possible  simultaneous  combination 
of  estimates  for  g0,  crx,  and  cry,  and  this  is  the  most  important  point.  These  small  errors 
probably  arise  in  part  because  the  equation  (5)  are  not  exactly  true. 

Table  3  shows  another  example  of  a  different  realization  and  alternative  input  values.  In  this 
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case,  E  passes  through  zero  at  ^  =  10.0  (highlighted  in  italics)  when  the  corresponding  ay 
and  9o  values  are  20.72  and  2.009  respectively.  These  estimates  are  again  very  close  to  our 
original  input  values  =  10.0,  a,  =  20.0,  and  9o  =  2.0.  Tables  2  and  3  clearly  demonstrate 
that  our  algorithm  determines  not  only  the  optimum  slope,  but  also  produces  estimates  of 
the  intrinsic  errors  associated  with  each  data  set.  Of  course  in  a  realistic  application  of  our 
algorithm,  we  do  not  need  to  produce  tables  like  those  just  shown  for  every  situation.  It  is 

adequate  to  let  the  computer  search  for  the  zero  crossing  point  of  E,  and  this  is  what  we 
normally  do. 

We  have  repeated  this  procedure  for  many  different  combinations  of  ax,  ay,  and  g0,  and 
for  many  different  realizations.  Figure  4  shows  a  collective  graph  of  the  ax/Ex  and  dy/£ 
values  which  we  have  estimated  compared  to  the  original  input  values  for  many  different 
realizations.  The  relationship  is  clearly  quite  linear  and  confirms  that  our  algorithm  is 
making  sensible  predictions. 

Figure  5  shows  differences  between  our  estimate  of  g0  relative  to  g0  after  subtraction  of 
unity.  We  clearly  produce  accuracy  of  ~  0.1%  to  1%  for  the  slope  of  the  line.  Thus  we  have 
demonstrated  the  validity  of  our  technique,  at  least  for  the  case  in  which  the  raw  data  and 
the  errors  are  normally  distributed. 

We  have  also  repeated  our  analyses  for  different  numbers  of  points  within  any  one  data  set. 
We  have  found  that  provided  there  are  greater  than  50  points  per  data  set,  we  produce 
error  estimators  for  ax  and  ay  which  have  absolute  errors  (Sdx  and  8dy)  of  less  than  about 
2.5%  of  E*  and  respectively  (see  Figure  4;  the  spread  in  dx/tx  and  dy/ty  is  about 
0.05).  However,  when  lesser  numbers  of  points  are  used  per  realization,  the  performance 
degrades  considerably.  Data  sets  of  20  points  become  substantially  less  reliable  with  regard 
to  estimating  ax  and  ay.  We  recommend  that  this  procedure  should  be  restricted  to  cases 
in  which  there  are  at  least  30  points,  and  preferably  50  or  more. 

We  should  finish  with  a  final  comment  about  multiple  minima.  As  noted,  we  have  created 
tables  like  Table  2  for  many  combinations  of  g0 ,  ax,  and  cry  (although,  as  we  have  said,  in 
any  practical  situation  we  allow  the  computer  program  to  find  the  zero  points  of  E).  We 
have  noted  that  for  some  cases  there  are  multiple  zero-crossings,  especially  for  cases  of  large 
<jy.  This  problem  is  worse  if  ax  >  cry,  so  we  recommend  that  if  users  have  any  a-priori 
knowledge  about  their  data,  they  should  choose  the  {a:;}  values  as  those  with  the  smaller 
ax.  Determinations  should  be  stopped  when  dx  exceeds  dy,  since  this  very  often  removes  any 
other  zero-crossing  cases.  If  no  a-priori  knowledge  is  available,  both  combinations  should  be 
attempted  (i.e.,  first  use  data  set  1  as  {x,},  and  then  use  data  set  2  as  {£;}). 
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ay/Iy  (calculated) 


a,/ 1,  ,  ay/Z,  (model) 


Figure  4:  Comparison  between  original  input  errors  and  estimated  values  obtained  from  our 
algorithm,  for  many  different  realizations. 
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Figure  5:  Differe: 
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1  1 
1  1 
1  1 
1  1 
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1  1 
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2  2 
2  2 
2  2 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
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4  4 

3  3 

2  3 

2  2 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
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40 

5  5 
3  3 
3  3 
2  2 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1  1 
1 


Table  4:  Number  of  zero-crossings  of  E  for  a  particular  combination  of  ax  and  cry,  for  £x  =  40 
and  =40. 


For  completeness,  we  show  in  Table  4  an  example  of  the  number  of  zero- crossings  of  E  for 
particular  combinations  of  crx  and  ay.  Clearly,  for  a  large  number  of  combinations  there  is 
only  one  minimum.  However,  there  are  cases,  especially  for  large  <ry  and  small  ax,  (where 
we  emphasize  that  these  are  the  “true”  intrinsic  errors),  for  which  there  are  several  zero- 
crossings.  As  already  noted,  we  recommend  dealing  with  cases  where  the  user  knows  that 
&x  <  o-y,  and  we  would  not  recommend  using  our  technique  if  the  user  detects  multiple  zero- 
crossings  in  the  region  dxfEx  <  dy/Hy.  In  fact,  most  of  the  “extra”  zero  crossings  associated 
with  the  cases  of  small  ax  and  large  ay  just  discussed  (see  Table  4)  arise  in  fact  when  our 
assumed  dx  values  exceed  dy.  This  is  a  nonsense  solution  and  is  avoided  if  we  do  not  carry 
our  iterations  in  Table  1  past  the  case  dx  —  <fy. 

Nevertheless,  There  are  clearly  a  wide  range  of  cases  in  which  we  can  uniquely  define  esti¬ 
mates  of  g0,  crx ,  and  ay.  If  the  user  can  determine  that  the  error  in  one  coordinate  is  less 
than  the  other,  and  uses  this  as  the  x  variable  (abscissa),  then  the  multiple  crossings  are 
rarely  a  concern.  The  first  zero  encountered  is  generally  the  correct  combination  of  g0 ,  dx, 
and  dy. 
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4  Conclusions 


A  new  method  for  least-squares  fitting  of  straight-line  data  has  been  presented  which  per¬ 
mits  optimal  fitting  and  simultaneous  determination  of  intrinsic  errors  for  both  the  abscissa 
and  the  ordinate.  It  assumes  that  the  original  data  and  the  intrinsic  errors  are  normally 
distributed.  The  algorithm  works  best  if  the  user  knows  which  of  the  two  variables  has  the 
smallest  intrinsic  error,  but  is  not  limited  to  this  condition.  There  are  occasions  when  the 
procedure  has  multiple  solutions,  but  these  are  relatively  rare  and  can  be  dealt  with  if  the 
user  exercises  due  care.  We  anticipate  that  this  process  should  find  application  in  a  wide 
range  of  situations. 
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